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Abstract 



We use Monte Carlo techniques and analytical methods to study the phase 
diagram of the M-component Widom-Rowlinson model on the bcc-lattice: 
there are M species all with the same fugacity z and a nearest neighbor hard 
core exclusion between unlike particles. Simulations show that for M > 3 
there is a "crystal phase" for z lying between z c (M) and Zd(M) while for 
z > Zd(M) there are M demixed phases each consisting mostly of one species. 
For M = 2 there is a direct second order transition from the gas phase to the 
demixed phase while for M > 3 the transition at Zd(M) appears to be first 
order putting it in the Potts model universality class. For M large, Pirogov- 
Sinai theory gives z<i(M) ~ M — 2 + 2/(3M 2 ) + .... In the crystal phase the 
particles preferentially occupy one of the sublattices, independent of species, 
i.e. spatial symmetry but not particle symmetry is broken. For M — > oo this 
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transition approaches that of the one component hard cube gas with fugacity 
y = zM. We find by direct simulations of such a system a transition at 
y c ~ 0.71 which is consistent with the simulation z c (M) for large M. This 
transition appears to be always of the Ising type. 
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I. INTRODUCTION 



The Widom- Rowlinson (WR) model, introduced in 1970 H as an ingeniously simple 
model for the study of phase transitions in continuum fluids (for an overview see Ref. |2j), 
continues to be, like its authors, a rich source of insights and analytical results in many 
(sometimes quite unexpected) areas of statistical mechanics. In this paper, dedicated 
with great pleasure to Ben Widom on the occasion of his seventieth birthday, we continue 
our study of a variation of the original model from two to M components on a lattice: hard 
core exclusion between particles of different species on nearest neighbor sites. 

This model was first considered by Runnels and Lebowitz who proved that when 
the number of components M is larger than some minimum Mq then the transition from 
the gas phase at small values of z to the demixed phase at large values of z does not take 
place directly. Instead there is, at intermediate values of z, z c < z < z d , an ordered phase in 
which one of the sublattices (even or odd) is preferentially occupied, i.e. there is a crystalline 
(antiferromagnetically ordered) phase in which the average particle density on the even and 
odd sublattices, p e and p Q are unequal. The average density, p(I), of species I — 1, . . . , M, 
on each sublattice is the same for each /, with p e (I) = M _1 p e and p {I) = M~ 1 p Q . The 
nature of the symmetry breaking is thus very different from that in the demixed phase at 
z > Zd where p e = p a = p but there exists one species, say I', for which p(I') > M~ x p. The 
origin of this crystalline phase is purely entropic. For z fixed and M large "it pays" for the 
system "entropy wise" to occupy just one sublattice without any constraint; since there are 
no interactions between particles on the same sublattice there are M independent choices 
at each site if we keep one of the sublattices empty. This more than compensates, at some 
M > Mo, for the "loss" of "fugacity energy" occasioned by keeping down the density in one 
of the sublattices. 

A natural question now arises, just how big does Mq have to be to see this ordered 
phase for M > Mq. It was shown in that on the square lattice Mq < 27 6 ; a ridiculously 
large upper bound. On the other hand a direct computation on the Bethe lattice |8|,[| with 
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g-neighbors gives M = [q/(q — 2)] 2 , which would suggest M ~ 4 for the square lattice, 
Mq ~ 3 for the cubic and M ~ 2 for the bcc lattice. Now it can be shown, using FKG 



inequalities, that Mq > 3 on any bipartite lattice ||10|| , but beyond that we have no simple 
or convincing argument for finding Mq. We therefore turned to Monte Carlo simulations. 
This gave on the square lattice M = 7 |||9| which is only about twice as large as the Bethe 
lattice prediction. This wetted our appetite to try the bcc lattice where q = 8. To our 
surprise we find here, using Monte Carlo simulations, that Mq does indeed equal 3, on the 
bcc lattice. 

While we have no clue of how to find rigorously the actual value of Mq or of z c (M), it was 
argued in || that for a given z and M large enough, the typical occupancy pattern on the 
lattice (ignoring the label I of the particles) should be like that of a one component lattice 



gas with nearest neighbor hard core exclusion. For the latter system Dobrushin |TT| proved 
the existence of a crystalline state. This implies that z c (M) should behave for large M as 
y c /M, where y c is the critical fugacity at which the one component hard cube gas (occupation 
at a site j excludes occupation at all eight neighboring sites) crystalizes. The value of y c 
for the bcc lattice, obtained by Gaunt [O using series expansion methods, is 0.77 ± 0.05. 



Using MC we obtained y c = 0.71 ±0.01 which is roughly consistent with Gaunt's value. Our 
result also agrees well with the values of Mz c (M) for large M being approximately .72 ± 02 
for M = 50 and M = 100. This provides solid evidence for the existence of a reentry phase 
transition in the temperature- magnetic field plane phase diagram of an Ising spin system 



with nearest-neighbor anti-ferromagnetic interactions on the bcc lattice [[L3| . 

We also find, as in ||, that for large M, z^(M) for the crystal-demixed transition can 
be computed via Pirogov-Sinai theory [TJ| yielding, 



M = z d + 2 - 2/?>z 2 d + ■■■ (1) 

which matches up smoothly with our MC results, see Fig. 3. It is easy to show that there is a 
demixing transition for M > 2 PJ5JH], the existence of sharp interfaces between coexisting 
phases, for M = 2, on the cubic lattice at large fugacity z was proven in ||. 



We next present results of our simulations and refer the reader to |J for a more detailed 
description of the model and additional references. 

II. RESULTS 

Our MC simulations were carried out on a bcc lattice of size 2 x S 3 = 2 x 22 3 with 
periodic boundary conditions. On an initially empty lattice we deposit particles chosen at 
random from the M components at fugacity z respecting the exclusion of different species 
occupying neighboring sites. We then sequentially update the lattice using a checkerboard 
algorithm resulting in a good vectorization. An update of a lattice site (4 )4 S \4 ) on one 
of the two simple cubic sublattices s, making up the bcc lattice, which is occupied by a 
particle of type I (I = indicating an empty site) is done as follows: We randomly choose 
a new trial particle of type I tr , where I tr can have any integer value between and M 
with equal probability. I tr = refers to an attempted removal of a particle I ^ from 
the lattice site, which is successful, if a number X randomly chosen with equal probability 
between and 1 is smaller than the inverse fugacity 1/z. When this occurs / gets the value 
0, otherwise it remains unchanged, I tr ^ refers to an attempted deposition of a particle 
of type I tr . If / = then it is successful if each of the four nearest neighbor sites is either 
empty or occupied by a particle of the same type (I tr ) and X < z. In this case / gets 
the value I tr , otherwise it remains unchanged. A direct replacement attempt of a particle 
/ 7^ surrounded by eight empty nearest neighbor sites is always successful. Typically in 
a simulation run after an equilibration of 2 x 10 5 Monte Carlo steps (MCS) we update the 
lattice 10 5 times, the configuration of every tenth step is taken for the evaluation of the 
averages. A typical run with 10 5 MCS took about 1 CPU hours on a CRAY-YMP. 

Let m(ii\ 4 , i^) denote the occupancy of a site, m(4 , 4 S \4^) = if the site 
(4. , 4 > 4 ) i s empty and m(i[ s \ 4 , 4 s "*) = 1 otherwise. As observables we took histograms 
Pl{4>c) of the order parameter C for the crystal structure and PL{4>d) of the order parameter 
4>d for the demixed phase in subsystems of size 2 x L 3 , 



5 



2L 3 

s=1 4 s) 4 s) 4 s) =i 



and 

1 



L£ jr f2m(4 s) ,4 s) ,4 s) )-ll (-1) 2 — 1 (2) 



3d 



■Max/JViCT) - p/M (3) 
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where Nl(I) denotes the number of particles of type J in a subsystem of size 2 x L 3 and p 
is the average overall density. 

A. Gas Crystal Phase Transitions 



For a given M the transition activity z c is found by finite size scaling techniques fll5| . |T6 
In particular the k — th moments of the order parameter distribution -Pl(0 c ), 

< <\>\ > L := [ <f> k c P L (<f> c )d<j> c (4) 



can be evaluated in subsystems of size 2 x L 3 , and from them the fourth order cumulant \W 
Ul , 

In a one phase region far away from a critical point the subsystem size typically can be 
chosen larger than the correlation length £, L » £ and the order parameter distribution is 
to a good approximation a Gaussian centered around 0, resulting in Ul — ► for L — > oo. In 
the two phase coexistence region far away from a critical point we can again assume L » £ 
and the order parameter distribution is bimodal resulting in Ul — * 2/3 for L — ► oo. Near 
the critical point however we have L << £, and using scaling arguments [16| the cumulant 
is a function of resulting for £ — > oo in the same value of C/# for all different L. This 
method allows the efficient determination of critical points by analyzing the cumulants for 
different values of z on different length scales L. Applied to our model we should see, for 
low values of z, when the system is in the disordered one phase region, Ul< > Ul for V < L. 
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Increasing z, we obtain, for large enough M, a crystal phase with TJu < Ul for V < L. 
Near z c we expect Uy ~ Ul for L' ^ L. This method for locating the transition fugacities 
was used in our previous studies |||| on the square lattice as well. 

For M = 3,4, 5, 10 we obtain in this way values of z c = 0.525±0.025, 0.3, 0.21, 0.085±0.01 
respectively, for M = 15, 17, 19, z c = 0.055, 0.047, 0.04 ± 0.001 and for M = 1000, z c = 
0.000725 ± 0.000025. Cumulants for the gas-crystal transition of the hard diamond and 
the M = 3-component system are shown in Figs. 1 and 2, respectively. The transition 
points are presented in Fig. 3 together with the asymptotic expression Mz c = 0.71. This 
corresponds to the value of y c on the bcc lattice y c = 0.71 ± 0.01 which we obtained using 
MC techniques. We also show in Fig. 3 a purely empirical fit of the l/M 2 -corrections for z c , 
z c = 0.71/ M + C/M 2 , to the MC data, with a value of C = 2 ± 0.5. The minimum number 
of components required for the existence of the crystal phase on the bcc lattice is M = 3. In 
Fig. 4 we show the approach of the critical fugacity Mz c to the limit of the hard diamond 
system. 



B. The Demixing Transition 

For M = 2 we obtain a direct transition from the gas phase to the demixed phase at a 
transition fugacity of z — 0.55 ± 0.02; the transition is second order with Ising exponents as 



on a simple cubic lattice studied by Dickman et al. |fL7|1 . 

For M > 2 we observe a direct first order transition from the crystal to the demixed 
phase. This transition was analyzed by studying the order parameter distribution Pl{4><i)- 
In the simulations we find a hysteresis region around Zd, going approximately between two 
values of z, say z\ and Z2, when increasing and decreasing the fugacity. In cases of a small 
hysteresis region with extent of less than \z% — z-i\ < 0.1, ^ was taken as the average of 
this region, Zd = (z\ + z^jl. In cases of a larger hysteresis region we located the transition 
fugacity Zd by finding the relative stability of one of the two phases during a simulation 
starting from configurations with both phases present in parallel slices extending over the 
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length of the simulation box, as described in ||. The resulting phase transition values of 
Zd are shown in Fig. 1. We note that with increasing number of components, the transition 
fugacities approach the exact asymptotic line M = z + 2 — 2/3z 2 + • ■ ■. 
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FIGURES 

FIG. 1. Cumulants versus z for the gas-crystal transition for M = 3. Symbols indicate 
different subsystem sizes, lines are for visual help. 

FIG. 2. Cumulants versus z for the gas-crystal transition for the hard diamond system on the 
bcc lattice. Symbols indicate different subsystem sizes, lines are for visual help. 

FIG. 3. Phase diagram in the M-z plane for a bcc lattice (MC). Full lines: Asymptotic lines 
for the phase transitions in the high fugacity region, M = z + 2 — 2/3z 2 , and for the transition in 
the low fugacity region, Mz = 0.71 ± 0.01, dashed line: transition line in the low fugacity region 
by fit of the 1/M 2 corrections for z c to the MC data, z c = 0.71/M + C/M 2 , with C = 2 ± 0.5. 
Symbols for MC: Transition points from the gas phase to the crystal phase (circles), from the gas 
to the demixed phase (squares) and from the crystal to the demixed phase (triangles). 

FIG. 4. Critical fugacities Mz c versus 1/M for the gas-crystal transition. The open symbols 
are results of the present Monte Carlo simulations, The black symbol (with error bars) on the 
y-axis indicates the critical fugacity of the hard diamond system from series expansions []l2|, the 
open symbol on the y-axis the critical fugacity from our MC simulations. Lines are for visual help. 
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Figure 1 (Nielaba and Lebowitz) 
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Figure 2 (Nielaba and Lebowitz) 
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Figure 3 (Nielaba and Lebowitz) 
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Figure 4 (Nielaba and Lebowitz) 



